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Abstract: One of the goals of power systems design and operation is to restrict the possible states of 
critical subsystems, during steady operation and during transients, to remain inside a certain bounded set 
of admissible states. Also, during transients, certain restrictions must be imposed on the time scale of 
evolution of the critical subsystem’s state. So, it is of interest to study if the state of the power system after 
a perturbation will remain bounded within some specified set of state variables and rates of change. If 
initially the system is in a steady state, it is of interest to study its behavior: it the power system will 
return to the original steady state or it will move farther away approaching to a new steady-state or to 
another attractor set, such as an oscillating pattern, or even suffer such a severe runaway that in some 
instant of time its physical integrity will be lost. 

The purpose of this paper is to apply certain asymptotic methods of nonlinear modal analysis (NMA), 
mainly but not exclusively to derive analytical formulae for thresholds of instability in mathematical 
models of power systems. 

After an assessment of some background aspects of NMA, certain modeling issues are discussed in to 
construct hybrid distributed parameters-lumped parameters mathematical description of power systems. 
In a first appendix to the main part of the present article, a comparison between NMA and Bunov- 
Galerkin methods is developed. 

As an example of hybrid modeling of an electrical-mechanical system, a lumped parameters model of a 
prime mover, a distributed parameters model of a connecting shaft and a lumped parameter model of 
loaded electric generator are combined. In a second appendix, from the hybrid model, a lumped 
parameters model is derived allowing an analytical approach to some details related with the interaction 
of the synchronous generator with power lines and loads. 

Two cases are considered with more detail in the main text. A static thermal instability of the “runaway” 
type in a nuclear reactor is analyzed by means of NMA methods and a global bifurcation is found. 

Then, nonlinear mechanical vibrations in a stretched shaft of uniform cross section are considered from 
the NMA point of view. Some advantages and limitations of the NMA methods are discussed. 


Keywords: power systems, nuclear reactors, turbine-shaft-synchronic electric generator assemblies, 
mechanical vibrations, nonlinear analysis, modal analysis, mathematical models, asymptotic 
analysis, instability thresholds, hybrid lumped parameters-distributed parameters mathematical 
models 


(A)-Introduction ! 

Power systems stability can be defined as a desirable property of any power system such that it 
is able to return to previous operating steady state after being subjected certain physical 
disturbances, with system variables bounded during transients, so the entire system remains 
intact (Kundur et al, 2004). 


' This article is a significantly expanded version of the original work prepared for the Queretaro 
International Conference on Dynamics, Instrumentation and Control of 2006. 
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In the present paper “power system” is used in a wide sense, as a general designation of 
electrical, mechanical, thermomechanical, electromechanical, hydraulic, fluid-thermo-chemical, 
systems where energy conversion and transport are produced in significant quantities, mainly at 
an industrial scale. 

One often used approach to study power systems stability is to construct a suitable mathematical 
model. Depending on the problem, this model is formulated as a system of ordinary or partial 
differential equations. To take into consideration time delays, sometimes the mathematical 
model is formulated using delay-differential equations of integral-differential equations. 

Thus, a dynamical system is obtained in a suitable state space. From a mathematical point of 
view, power systems are highly nonlinear dynamical systems. Besides, the mathematical model 
must consider that these systems operate in a constantly changing environment: power inputs 
and outputs as well as important operating parameters change continually, and the stability 
behavior of the system depends both on the initial operating conditions and the nature of the 
disturbances. 

When the state space is of finite dimension, the approach is by lumped parameter systems. 
When fields, functions of time and space, are used in the mathematical model, the approach is 
by distributed parameter systems. 

The best-known mathematical method to study stability problems of lumped parameters 
dynamical systems, given as systems of ordinary differential equations, is the First Method of 
Liapunov, also known as local linearization, or stability in the presence of infinitesimal 
perturbations. It allows us to study strictly only the stability of the steady state of a linear 
version of the system of nonlinear ordinary differential equations mathematical that describe a 
dynamic system. However, with the aid of a theorem due to Hartman and Grobman, its scope 
can be extended to study the stability of non-linear systems in a small enough neighborhood of 
the steady state. 

Although the First Method of Liapunov can be generalized to mathematical models of 
distributed parameter systems including systems of partial differential equations, it can't be 
applied to stability problems in the large (that is, the response to non-infinitesimal 
perturbations), like the ones posed by certain mathematical models used in the design, operation 
and control of power systems. 

The analytical methods of non-linear modal analysis (NMA) allow us to tackle some of these 
problems, both in lumped and in distributed parameter systems. They can be considered an 
extension of Liapunov's First Method to cope directly with the non-linear terms in the dynamics. 
These methods complement both, numerical methods ab-initio, and the results obtained after 
applying some summarizing function criterion to determine a region of stability, like Liapunov's 
functions (Liapunov’s Second or Direct Method) and the like, that are often applied in 
theoretical discussions of power systems dynamics and control. 

The NMA analytical methods were developed, by W. Eckhaus and others, to study finite 
amplitude instability problems in fluid mechanics. During the sixties and the seventies, the 
NMA methods were applied to the study of stability problems in distributed parameter systems, 
including mass transport processes and chemical reactions in continuous flow industrial 
chemical reactors (fluid-thermo-chemical systems) and other systems of interest from an 
engineering standpoint. But NMA methods for distributed parameter systems do not seem to 
have attracted the interest of people working with mathematical aspects of other kind of power 
systems. It was Liapunov’s Direct Method that attracted their attention. This was mainly due to 
the capability to analyze the finite amplitude behavior of nonlinear equations without ever 
solving them. Almost all the work was done beginning directly with lumped parameters models 
of the power system. 
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In modal analysis the different fields that give the space-time evolution of a distributed 
parameters system with a bounded space domain, are expanded in series of eigenfunctions of a 
certain suitably chosen linear operator, including the boundary conditions of the problem. These 
eigenfunctions (spatial modes) depend only on the space coordinates and are defined in the 
bounded domain. In each expansion corresponding to a given field, these spatial modes are 
weighted by unknown time dependent modal amplitudes. The modal amplitudes may be 
determined so that the series expansions represent a solution of the dynamic field equations 
verifying the initial and boundary conditions in the bounded space domain. The original 
nonlinear partial differential equations with time and space as independent variables are 
substituted by an equivalent system of nonlinear ordinary differential equations for the unknown 
mode amplitudes, with time as the only independent variable. From the initial conditions 
verified by the fields, a set of initial conditions for the mode amplitudes are obtained. Thus, a 
complex space- time field dynamics is reduced to the study of the evolution of a representative 
point in the space of mode amplitudes. Often it is possible to work with a relatively small 
number of mode amplitudes, after reducing the original high order dynamics to a low order one, 
applying the so-called inertial manifold theory (Malinietski, 2005). 


The purpose of this paper is to apply certain asymptotic methods of NMA, mainly but not 
exclusively to derive thresholds of instability in distributed parameters power systems. 

After an assessment of some background aspects of NMA, some modeling issues are discussed 
in order to construct a hybrid distributed parameters-lumped parameters mathematical 
description of the power system. 

In a first appendix to the paper, a comparison between NMA and Bunov-Galerkin methods is 
developed. 

A static thermal instability of the “runaway” type for a nuclear reactor is studied by NMA 
methods. 

Next, nonlinear mechanical vibrations in a stretched shaft of uniform cross section are 
considered from the NMA standpoint. 

In a second appendix, from the hybrid model description of the power system, a lumped 
parameters model is derived allowing an analytical approach to some details related with the 
interaction of the synchronous generator with power lines and load. This allows an easier 
analysis of some aspects of sub-synchronous resonance phenomena. 


Approximate analytical formulae that relate the instability thresholds and oscillation periods 
with the main parameters of the mathematical model of the system are derived. 

Some advantages as well as some limitations of the NMA methods are highlighted in the final 
discussion and conclusions. ” 


(B)-Asymptotic methods of Nonlinear Modal Analysis 


? Today it is often easier to make an ab-initio numerical analysis of transient dynamics and stability of 
power systems, even without the introduction of the mode concept. This is the contemporary trend, 
fostered by the increase in computer power and the development of high-quality numerical algorithms 
(Conca and Gatica, 1997, De Silva, 2005). However, only through digital simulations it could be difficult 
to show clearly enough the manifold of effects, on the studied dynamics, of changing the physical 
parameters of the system. So, a simplified analytical approach may be of interest, both as a guide for the 
design of ab-initio digital simulations, and for the interpretation of experimental results. 
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(B.1)-Nonlinear modal analysis in a generalized sense 


In a generalized sense modal analysis of partial differential equations is a procedure that allows 

(Attlee-Jackson, 1991; Zauderer, 2006): 

(a) The expansion of the solutions (fields) of partial differential equations in time and space 
variables as a series of given space functions weighted by time functions known as modal 
amplitudes. 

(b) The derivation of evolution equations for the mode amplitudes (a set of ordinary 
differential equations). 


As the space functions are known, by means of modal analysis a complex space- time field 
dynamics is reduced to the study of the evolution of a representative point in the space of mode 
amplitudes. 

If the original partial differential equation is nonlinear, the ordinary equations for mode 
amplitudes are nonlinear also: in this case we have nonlinear modal analysis. 

An often-used procedure to obtain the equations for mode amplitudes is the application of 
Galerkin’s method (Gershenfeld, 1999). In general, in the evolution equations thus obtained, 
modal amplitudes appear coupled, even in the linear approximation. 


(B.2) Nonlinear modal analysis in a restricted sense 


In a restricted sense, modal analysis is a procedure applied to the different fields that give the 
space-time evolution of a distributed parameters system with a bounded space domain D, 
expanding the fields in series of eigenfunctions of a certain suitably chosen linear operator L| |. 


The domain U of this operator, which is a set of space dependent functions defined in D, is so 
defined to include the boundary conditions of the problem. These boundary conditions must be 
time independent, although they can be non-homogeneous. 

The abovementioned eigenfunctions (spatial modes) depend only on the space coordinates and, 
as already said, are defined in the bounded domain D. 

In each expansion corresponding to a given field, these spatial modes are weighted by unknown 
time dependent modal amplitudes. 

The modal amplitudes may be determined so that the series expansions represent a solution of 
the dynamic field equations verifying the initial and boundary conditions in the bounded space 
domain. 

The original nonlinear partial differential equations, with time and space as independent 
variables, are substituted by an equivalent system of nonlinear ordinary differential equations 
for the unknown mode amplitudes, with time as the only independent variable. One advantage 
of this restricted kind of modal analysis, compared with the generalized one, is that in the 
obtained evolution equations the modal amplitudes always appear decoupled to the linear order 
(Suarez Antola, 1997; 2005 a; 2005 b). 

From the initial conditions verified by the fields, a set of initial conditions for the mode 
amplitudes are obtained. 

Thus, as restricted modal analysis is part of generalized modal analysis, also in this case a 
complex space- time field dynamics is reduced to the study of the evolution of a representative 
point in the space of mode amplitudes. 

Often, it is possible to work with a relatively small number of mode amplitudes, after reducing 
the original high order dynamics to a low order one: this is the subject of inertial manifold 
theory. 


International Conference on Dynamics, Instrumentation and Control, Querétaro, Mexico, August 2006. 


A necessary condition to be able to work, at least asymptotically, with a small number of time 
dependent mode amplitudes is that the eigenvalues of the considered linear operator be widely 
spaced. 


(B.3)-Relation between nonlinear modal analysis and Liapunov’s First Method. 


Now, there is a close connection between the classical Liapunov’s First Method to study 
stability problems and nonlinear modal analysis. 

The First Method of Liapunov, also known as local linearization, or stability in the presence of 
infinitesimal perturbations, allows us to study strictly only the stability of the linear version of 
the dynamic system. However, with the aid of the well-known theorem due to Hartman and 
Grobman (Saaty and Bram, 1964; Attlee-Jackson, 1989), its scope can be extended to study the 
stability of non-linear systems in a small enough neighborhood of the steady state. But it can't 
cope with stability problems in the large (that is, the response to non-infinitesimal 
perturbations). However, the design, operation and control of engineering systems often pose 
stability problems in the large. The analytical methods of nonlinear modal analysis (NMA) 
allow us to tackle some of these problems. 

So, seen from this standpoint, NMA can be considered an extension of Liapunov's First Method 
to cope directly with the non-linear terms in the dynamics. These methods complement both, 
numerical methods ab-initio, and the results obtained after applying some summarizing function 
criterion to determine a region of stability, like Liapunov's functions (Liapunov’s Second or 
Direct Method) and the like (Saaty and Bram, 1964), that are often applied in theoretical 
discussions of engineering control problems. 


(B.4) Foundations of restricted modal analysis 


As an example, let us consider the following nonlinear evolution equation for the fieldu(t, x ) : 


yu} Nlul+f [B.4.1] 


The operator Llulis linear, and N [u is a nonlinear one. Both in the linear and in the nonlinear 
operator appear the field u(t,F ) and its space derivatives, multiplied by suitable coefficients 


that can be functions of time and space coordinates. 
We assume that in the nonlinear operator there are no linear terms (all linear terms appear in 


Lu] by construction). 
The term ff is a source or external forcing function f (t, r ) ; 


Furthermore, we suppose that the zero solution is our reference solution of equation [B.4.1], 
although in the original problem we may be faced with a steady state field with spatial 
variations. If this is the case, the equations are manipulated until an equation like [B.4.1] is 
obtained (this can be always done). 


Then we can expand the nonlinear operator in a Taylor series, using the Fréchet derivatives 
introduced in functional analysis (Griffel, 2002): 


Nu] =5-D*[Ouu}+ = D'osuuu}. [B.4.2] 


The operator D’(0;u,ulis bilinear, D*[0;u,u,u] is trilinear, and so on. 
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We assume that the linear operator has simple eigenvalues A, with eigenfunctions @, such that 


the set of eigenfunctions is complete in the space of functions that has the considered field, at 
each instant, as one of its elements. 


Now, let us consider the adjoint operator L’ [u | corresponding to Llu], its eigenvalues A, and 
its eigenfunctions @' p(X): aa "|= 1, °D" p [B.4.3] 
The set of eigenfunctions of the adjoint operator is also a complete set in the function space to 
which the field u (t, 7) belongs. 


Each eigenfuction 2, is orthogonal to the eigenfunctions @, of Llu] that corresponds to the 


other eigenvalues: (9, P,) =0 D#q [B.4.4] 
It is always possible to put: (9; Pp) =A] 
Next, we define the mode amplitude A w(t) by: 
A, (t)=(9;,.u) = | 9; ( 7)-w(F)-d7F =0 [B.4.5] 
D 


Here w(F ) is a positive weight that may be identical to 1. Then, with the modal 
amplitude A_, (t) given by equation [B.4.5]:  u (t, 7) = SA, (1) p. (7) [B.4.6] 


Making the internal product of both members of [B.4.1] with 2, we obtain: 


* 0 * * * 
(e; =) =(9},,Lu)) +(9;,Nlu) +(9;,f) [B.4.7] 
. u\ d;. \_ 4A,() 
But: (0; ) = = \Pr#) = oat [B.4.8] 
(9, Lu)) =(L lo; fu) =1,-A,{(t) [B.4.9] 
; iD Apps ioe ee 
(y;,.[u}) = (9, ,D* [0;v,w]) + = (g;,,.D°[0;u,u,u)) oe [B.4.10] 
Now we substitute the expansion [B.4.6] in equation [B.4.10], and after several calculations it 
follows: 
(g,, Nu) = ae -A,-A, + YOr.2 -A,-A,-A, +... [B.4.11] 
qg.r=l q.r,s=l 
Here, by definition: Pia = (9, ,D* lo: 2, 9) [B.4.12] 
Oars = (9;,.D*[0: 2.9) [B.4.13] 
Some of the symbols @,, |, @,.4.,,> «++ a8 Well as the eigenvalues 2, may be functions of time. 


Putting together equations [B.4.7] to [B.4.13] and defining the projection of the source term 
onto p,(F )by ts (*)=(9%, vA Ns we derive a numerable infinite set of ordinary differential 


equations in the mode amplitudes: 


dA, E 
= A, As o Dens 


-A,-A,+ ) 6 
dt q.r=l a 24 


P»Gst 8 


-A,-A,-A,+..+f,(t) p=12.... [B.4.14] 
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These are the dynamical equations for the so called restricted nonlinear modal analysis, 
restricted because in the linear approximation all the modal amplitudes appear uncoupled: 


dA 
sed =A, -A,(t)+ f,(t) DP = 12525 


As will be shown below, this approach to NMA can be applied to equations formulated in terms 


rag 


of Yaad ) and with cross derivatives of time and space variables, as well as to evolution 
t 


equations formulated in terms of several interrelated fields (state variables of the distributed 


parameters system) u, (t, 7) JU, (t, ae dat, (t, f) 


In the Appendix (F.A) a summary of the methods weighted residuals is given. The methods of 
weighted residuals are used to convert a partial differential equation with time and space 
position as independent variables, into a system of ordinary in the time variable only. 
Bubnov-Galerkin method is a particular case of these weighted residuals methods. 
o ae ae ; 
Working with the nonlinear wave equation — u(t, x) =c? ——u(t,x)+¢| ———u(t,x) | , it is 
r) rt? ( ) Ox? ( ) S r) t 6x? ( ) 
shown there that Bubnov-Galerkin method is closely related with the procedure developed here 


as a foundation to restricted NMA. 


(B.5)-The meaning of mode coupling and the advantages of restricted modal analysis to 
study stability problems. 


As the superposition principle does not apply in the case of nonlinearly coupled mode 
amplitudes, one may ask about the utility of nonlinear modal analysis. 

The nonlinear interaction consists in the exchange of energy between different modes. 

Thus, even if the initial field may be represented by only one or a few modes, the nonlinear 
processes tend to activate other modes as time progresses. 

In a Fourier setting this would be a generation of harmonics or sub-harmonics. 

Besides the possibility of creation of spatial patterns of smaller or larger wavelengths, there are 
other important effects: self-damping of a growing mode and inter-mode suppression. 


The self-damping may appear in a mode that is unstable according to the linear theory 
(A, >). If this mode is the only one that is perturbed from its steady state its amplitude first 


begins to grow exponentially: A, (=A, (0)-exp(Z, -t) When, ,, <0, the quadratic 


term tends to slow down the growth of the amplitude (self-damping). 


If the initial perturbation affects a certain number of different modes and if initially many of 
these grow according with the linear approximation, the growth of one mode could tend to 
damp the growth of the other modes (inter-mode suppression), either because it is the fastest 
growing mode or because its initial amplitude is much larger than the initial amplitudes of the 
other modes. 
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This occurs when there is some competition between modes, for example a competition for a 
given amount of energy. 


Both self-damping and inter-mode suppression allow a partition of the space of mode 
amplitudes into basins of attraction of equilibrium points, leading to bi-stability, run-away and 
other interesting dynamic behaviors. 


When there are dominant eigenvalues, the evolution equations cast according to the restricted 
procedure of nonlinear analysis are particularly suited to: 

(a)-Study mode slaving in Haken’s sense (Haken, 1983). 

(b)-Apply Eckhaus’s method (Eckhaus, 1965) to determine sufficient conditions for instability. 
Analytical formulae can be derived to characterize the dynamic behavior of a small number of 
dominant modes (applying asymptotic methods, see Haken, 1983; Nicolis, 1995). The 
understanding thus obtained serves as a guide to the design of in-depth digital simulations, using 
modal equations or directly by discretization of the original nonlinear field equations. 


(C)-Hybrid distributed parameters-lumped parameters mathematical description 
of power systems 


The definition of power systems stability can be applied either to an interconnected power 
system as a whole or to one of its components. For example, in a complex electromechanical 
system, with many generators of electric power, transmission lines and electric motors, one 
generator or one motor may lose stability without producing a cascading instability in the main 
system (Kundur et al, 2004). 

In this case the loss of stability of the component can be studied considering the component as a 
system by itself and assuming the rest of the interconnected power system as its environment. 
When an analytical approach to a stability problem is desired, it is often necessary to work with 
a hybrid of distributed parameters with lumped parameters models. 


If a rotating machine is the focus of interest, from a mechanical engineering point of view it is 
stable if its rotor performs a pure rotational motion around an appropriate axis at a required 
rotational speed and its motion is not accompanied with other modes of vibration of the rotor, 
its component parts, or other stationary parts of the machine, or, if such vibrations take place, 
their amplitudes do not exceed admitted, acceptable values. A stable rotating machine is 
immune to external perturbing forces: any random perturbation cannot drastically change its 
behavior. Such perturbations cause only a transient decaying process leading to a previous 
regime of performance, or to a new one, which is included in the acceptable limits (Muszynska, 
1984). 


In the large turbine-generator rotor trains of nuclear power plants, lateral vibrations, torsional 
vibrations, and stator vibrations can be found. 

Lateral vibrations can be due to: classic (Laval) shaft mechanical unbalance, local heating that 
leads to a thermal unbalance equivalent to a mechanical unbalance, a bow of the generator rotor 
due to unsymmetrically distributed friction forces in circumference direction between the 
cupper windings and the generator steel rotor, different stiffness values of the turbo-generator 
rotor in two directions, effects on the steam flow (across the low pressure section of turbine) 
produced by variations in the inlet pressures of the turbine’s primary coolant circuit related to 
changes in the secondary cooling water temperature or in electrical loads. 
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Torsional vibrations can be due to due to electrical faults (synchronous or super synchronous 
vibrations), sub synchronous resonance and blade vibrations in last stage of the turbine due to 
fluid-structure interactions. 

Synchronous or super synchronous torsional vibrations due to due to electrical faults (for 
example two-phase and three-phase short circuits) are produced by oscillating air gap torques of 
the network frequency and twice the network frequency. As consequence, torsional vibrations 
with resulting internal torques and shear stresses in the shaft elements are excited. 

Sub synchronous resonances are oscillations that sometimes appear both in the electrical and in 
the mechanical systems of a power station. They are due to the interaction between the 
turbogenerator and its power line when it is compensated by series capacitors. An electrical 
resonant circuit can be triggered by a suitable fault in the electrical system. An exchange of 
energy can occur between the shaft and the inductive and capacitive elements of the electric 
circuit, generating low frequency (well below the synchronous frequency) electrical torques in 
the generator air gap that act on the generator rotor. If these frequencies are near one of the 
(lowest) natural torsional frequencies of the shaft line, strong resonant vibrations may appear in 
the shaft assembly. These vibrations are in turn transmitted into the electrical system by the 
electromechanical coupling. 

Blade vibrations can be produced either coupled with the torsional vibrations of the rotor system 
(in the last stage blades of the low-pressure stage of the turbine group) or excited by the steam 
in the turbines of the shaft train. 


In a turbo-generator node of a power system, some stability issues during the transfer of power 
between the turbine and the generator can be related with vibrations and delays produced by 
wave propagation in the shaft. 

Some flexible shafts linking a mechanical generator with a mechanical load, in the framework 
of a control system, sometimes show sustained parasitic torsion oscillations. When the sensor of 
a feedback system whose output modifies the actuator’s output, measures the local response of 
the shaft, the propagation time of the torsional waves from the actuator location to the 
measuring instrument location introduces a delay in the feedback loop. If this time delay is big 
enough, it can cause instability. 

A very coarse lumped parameters mathematical model allows us to see this. Let us suppose that 


the mechanical power P,, (1) to be delivered by the prime mover must be a certain constant 
power P*. 


Let us further suppose that, due to a simple feedback loop, P,, (t) verifies the following delay 


differential equation: —P(t)=P,"-P, (t-t,) [C.1] 


m m 


Here T. is a characteristic response time of the prime mover with feedback. 


m 


Introducing p(t) =P (=P, in [D.1] we have: T. < p(t) = —p(t =f) [C.2] 
t 


In this case (Hayes, 1950; Driver, 1977; Minorsky, 1983) if t,<T7. then lim p(t)=0: as 
t 


+00 


desired, the mechanical power P_ (t) tends to P,” as time goes on. But when the inequality is 


reversed, p(t) oscillates with increasing amplitude. 
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This means that the power P, (t) oscillates with increasing amplitude around P * as time goes 


on until the mathematical model cannot be applied any more due to the always present 
nonlinearities. 


The torsion, flexural or combined torsion-flexural waves in shafts can be studied by suitable 
(dissipative) wave equations (a distributed parameter model) with boundary conditions coming 
from a mechanical generator in one end (prime mover, for example a steam turbine) and a 
mechanical load (mainly of electromagnetic origin) in the other end. In this case both, prime- 
mover and mechanical load can be described by a lumped parameter model that details their 
interactions with the environment. 

This environment includes the rest of the power system: the thermal-mechanical subsystem 
(nuclear reactor or boiler, steam generator, fluid pipes, pumps, etc.) and the electrical power 
subsystem (synchronous generator and its transmission lines, transformers, interconnections 
with other generating nodes, electrical loads, etc.). 

When the mechanical effects of hot spots in the rotors cannot be neglected (shaft differential 
heating or Morton effect), the corresponding thermal equations and their coupling with rotor- 
dynamics must be specified. Besides, the mathematical model should consider the effects (often 
nonlinear) due to the stiffness and damping capacity of the fluid film bearings used in heavy 
rotating machines. 


The following schema (Figure 1) highlights the genesis of mechanical oscillations as 
consequence (and part) of the energy exchange between the mechanical and the electrical power 


systems. 
Oscilati 
Mechanical System 
sausuaumanaen — Electrical System 
Mechanical oscillation 
M, 


hu 


Fig. 1: Mechanical-electrical oscillating energy exchange (modified from a Siemens report) 


A simplified hybrid mathematical model of this exchange is the following generalization of the 
classic Stodola lumped parameters vibration model. * It will be built as an example of hybrid 
modeling, done with the intention of applying NMA methods. 


3 This model was originally posed to analyze the possibility of having the period of the alternating current 
in a synchronous generator near to the period of natural vibrations of the common shaft that connects the 
generator rotor with a turbine prime mover. Stodola (1924) studied a two degree of freedom lumped 
parameter model. The shaft connecting the turbine to the generator was considered as a massless linear 
torsion spring loaded with two inertia moments at its ends, one due to turbine and the other due to 
generator rotor. All mechanisms of energy dissipation were neglected in Stodola model, as well as the 
always present unbalance due to some degree of material heterogeneity or manufacturing errors. 
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Both the generator rotor and the turbine are considered from a lumped parameters point of view, 
being O, and ©, their respective inertia moments, being >, and ¢g, their respective rotation 


angles. A uniform shaft of circular cross section and length Lconnects these rotors: only 
torsional vibrations will be studied. 4 
The shaft is considered from a distributed parameters point of view, so its equation for torsional 
oscillations is supposed to be the following linear and dissipative wave equation: 
2 2 3 
Ee (eee eee Ce [C.3] 
Cc; Ot Ox Otox 


Here O(t, x)is the rotation angle of the shaft. If the shaft axis is given by the 


interval0<x<L,then: O(t,0)=¢(t)  [C4al O(t,L)=¢,(t)  [C.4b] 
3 


The dissipative term 7, 5 O(t,x) follows from the assumption that local shear stress S is 


Otox 


0 
related with local shear strain y by the constitutive relation: S=/7+T, a Y [C.5] 
t 


If G is the shear modulus and p is the density of the shaft material, the shear wave velocity c, 


G 
is given by: e =— [C.6] 
fe) 
If J , 18 the polar moment of the shaft, the local torsion moment is given by: 
Wer | 2 ae dee 0s) C7 
Ox Otox 


At the turbine-end of the shaft, if M ae )is the assumed known mechanical torque (prime 


mover torque) produced by the turbine blades: 


a 

©, Fras (t)=—-M_,,(t)+M, (1,0) [C.8] 

From [C.4 a], [C.7] and [C.8] we obtain the boundary conditions in the turbine end of the shaft: 
0 a on 

GJ, [2 otro) 45 ae O(t, 0)}- 0, De —, 9(t,0) =M.,,,(t) [C.9] 


At the generator-end of the shaft, ifM ont )is the electromagnetic torque acting on the 


generator rotor: 


0, 4, (t)=—M, (1,L)+M.,, (t) [C.10] 


From [C.4 b], [C.7] and [C.10] we obtain the boundary conditions in the generator end of the 
shaft: 


4) eo a 
61, Zot) a age 1)|+0, az Olt L)=M,,,(t) [C1] 


General initial conditions are: 0(0, x) =0, (x) [C.12 a] = 0(0.x) =O, (x) [C.12 b] 


4 Here, as it was done in the above mentioned Stodola’s model, the always present unbalance (due to 
some degree of material heterogeneity or due to manufacturing errors) will be neglected. 
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If the prime mover torque and the electromagnetic torque are known, equation [C.3] with non- 
homogeneous boundary conditions [C.9] and [C.11] and initial conditions [C.12] can in 
principle be solved by analytical (approximately, using asymptotic methods when the 
dissipative term is small enough) or numerical tools. 

If an analytical approach is intended, the often-small dissipation of torsion energy allows a 
perturbation analysis. 


Let us define the following dimensionless space € and time t variables: x=LE& t= <r T 
T 
In these new independent variables, equation [C.3] can be rewritten thus, where ¢ = or Fo ; 
2 2 3 
ee té g 59 [C.13] 
OT 0g Or0g 


The boundary conditions [C.9] and [C.11] can be recast as follows: 


é Oo 0, | a 2 eon 
[ Zo(e.0)6 2 0(6.0)-[ 52} Z0(.0)-| Jat) [C.14] 


a e° eo, \e Aa 
[Zo(e.)ee sgtet)}o{ 9} 0(6)=| Je (zr) [C.15] 


If ¢= “is small relative to 1, a perturbation analysis can be done, assuming M_,,, (z) and 
L 


M en (T) as known and considering the presence of the small dimensionless parameter in the 


derivative of highest order, both in the partial differential equation and in the boundary 

conditions. 

Suitable methods of perturbation theory are available when dissipation cannot be neglected 

(Holmes, 1995; Zwillinger, 1997; Zauderer, 2006). 

If dissipation can be neglected, in [C.13], [C.14] and [C.15] we put ¢=0 and obtain the wave 
2 2 

0 =p 0 . 

OT 0g 


equation @ together with the boundary conditions: 


Now, the method of Mindlin and Goldman (Meirovitch, 1967; Clark, 1972), given the initial 
conditions, allows us to find the solutions in this simplified case. An alternative approach is to 
construct a solution by the method of finite Fourier transforms (Zauderer, 2006). 


However, if a certain understanding of the exchange of energy between the mechanical and the 
electrical subsystems is the goal, some modeling of the turbine and the electric generator must 
be done. A minimum of modeling of the synchronous generator, its transmission lines, electrical 
loads, following well established procedures (Woodson and Melcher, 1968) is necessary to 
describe the energy exchange in the generator end of the shaft and the torsional oscillations 
related with this exchange, even assuming a constant torque due to the prime mover in the 
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turbine end of the shaft. This broadened modeling, depending of its details, introduces 
nonlinearities in the mathematical problem. 

From the hybrid distributed parameters-lumped parameters model built in this section a simple 
lumped parameters model can be derived and easily studied, as shown in the Appendix (F.B). 


(D)-Examples of NMA approach to static and dynamic instabilities in power 
systems 


Here, in (D.1) an example of a static instability of the “runaway” type for a nuclear reactor is 
given (this example is taken from Sudrez-Antola, 2005 a). 

Then, in (D.2) nonlinear mechanical vibrations in a shaft of uniform cross section are 
considered from the NMA point of view. 

Approximate analytical formulae that relate the instability thresholds and oscillation periods 
with the main parameters of the system are derived. 


(D.1)- An example of a static instability of the “runaway” type in a mathematical model of 
a nuclear reactor 


Let us consider a homogeneous bare nuclear reactor, whose extrapolated core fills a region B 
(for an introduction to nuclear reactors and nuclear engineering see Lamarsh and Baratta, 2001. 
For nuclear reactor analysis more in depth see Duderstadt and Hamilton, 1976 or Bell and 
Glasstone, 1971. For nuclear reactor physics see Lamarsh, 1966 or Stacey, 2001). 

This example is posed only to show as directly and as simply as possible how threshold 
amplitudes for instability can be derived from the methods of nonlinear modal analysis. Because 
of that it is not intended to be a realistic model to be applied to some practical problem. For that 
reason, we shall work with only one group of neutrons, and with only one group of delayed 
neutron precursors. We suppose that the neutron field flux, the precursor’s field and a feedback 
variable verify the equations of evolution: 


FA £ 
ct = M’ V7¢+([1- 8] (K, + kw kw’) —1) 6+ 4, +G [D.1. la] 
0b, _ Blk, +kw-k,w’) 
es ee oo: [D.1.1b] 
Ot L 
ma Ow _ Kb [D.1.1c] 
Ot 


Here t is time, r is position vector, d(t,F ) is neutron flux, 4 is a characteristic lifetime of 
neutrons, f is the delayed neutron fraction, @, (t,7 ) is a flux proportional to the density of 
delayed neutron precursor atoms, ty=//A is a characteristic time scale of precursor decay (with 
decay constant i), M * is the quotient between the diffusion coefficient and the macroscopic 
absorption cross-section, k is the multiplication factor and S,is a source term due to the 


presence of external neutrons that are being emitted with this emission rate, w is a feedback 
variable, 7, is a characteristic time scale of evolution of w, and the parameter K links the 
neutron flux with the rate of variation of w. 


Furthermore, we supposed that there are some feedback mechanisms that make the 
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multiplication factor k a polynomial function of w: 

k=k,+kw-k,w [D.1. 2] 
Herek,, k, and k, are positive parameters. As consequence, if w increases from zero, k first 
increases and then decreases. The boundary conditions are dlt, 7, ) =0 ,¢, (r Li ) =0, and 


w(t, i, ) =0 , forevery t and every 7, € OB (the boundary of region B). 


Now, our purpose will be to study the stability of the stationary solutions of equations [D.1.1a], 
[D.1.1b] and [D.1.1c]. 
Let us suppose that the source of external neutrons vanishes, so that g, (7 ) =0 P10 (7 ) =0, and 


Wo (F ) =(), is a possible steady-state solution (zero power solution). 


However, we may introduce a suitable distributed source of external neutrons, during a while, to 
produce a perturbation in the field variables relative to the strictly zero-power solution. After the 
desired perturbation is produced, the source will be removed. 


In this case the linear operator Ap [0] obtained fixing w at its steady state value w = w jis: 
M? a 


Adld]=—-V°b+(|- B)-4 


The nonlinear operator Ny [¢; w] now is given by N, [¢; w] = ( =5) k,w@ ( =) kw. 


Then equation [D.1.1a] may be cast as follows, in a form suitable to begin with NMA methods: 
ag > ‘ 1 1 
— = A.|d]+ N,[¢;w]+—-¢, +—S, 
Ot a ¢ 
In this case the eigenfunction- eigenvalue problem for the operator Ald] with homogeneous 
boundary conditions for ¢ is equivalent to the eigenfunction- eigenvalue problem for the 


operator — V* [d| with the same boundary conditions for@ . 


The eigenfunctions are the same for both operators, and the corresponding eigenvalues are 
related by a straight line. 
As consequence, to study the stability of the zero-power solution, we represent the three fields 


as series of eigenfunctions of the operator -v*| | with homogeneous Dirichlet’s boundary 


conditions in the boundary of the domain B : 


o(t.7)=Y 4,09, (*) (D.1.3 a] 


n=0 


w(47)=>B,0¢,(7) [D.1.3 b] 


n=0 


p, (7) =YC,04,(¥) [D.1.3 c] 


Here -V’ ¢, (7) = (7) , and the eigenvalues can be ordered 0 < up < Ly < Ly <u 


The eigenfunctions are orthogonal and can be normalized. 
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Then ( (¢ wP m) =I ¢ ,( =6,,, being 6,,,, Kroenecker’s delta. 


Substituting the ansatz [D.1.3] in equations [D.1.1] and projecting onto each eigenfunction @ ig 


the following infinite system of nonlinear ordinary differential equation is obtained: 


0A, =([1-A]k (14 2x8)) A, +f +[1-p]k, y Lynn AnBy -[1- B] ke y Lng 4yB,B, 
m,n=0 m,n,q=0 
C4) 
Ty 
[D.1.4 a] 
OB, 
ae ne [D.1.4 b] 
eC. ie kA, +k, ae ann B, —K, yi ring Am BB, |- Cy [D.1.4 c] 
“at t rele m,n,g= i : : 
p=0,1,2,3... 


Here Fyn, =[9,GnGAV + Tong =]%,GnGrG,4V and, =| $,(7)S o(t,F)dV 


If A. (0) =( for every mode index p=0,1,2,3,....and at least one of the forcing 
terms S$ a 0, then a perturbation will be produced from the zero-power solution. 


Once every projection S o (t) has vanished, the perturbation is already completed, and the next 


task is to determine its future evolution. 


Will the flux return to zero everywhere, or at least remain bounded and relatively near zero? 
Or will it suffer a severe runaway, either approaching a new steady-state solution or even 
reaching values that put in danger the physical integrity of the reactor? 


To continue, let us suppose that at zero-power the reactor is sub-critical. This means that the 
inequality k, <1+M°y; is verified, and because LL, >p, if p #0, all the coefficients of 
the linear terms are negative in equation [D.1.4 a]. 


Let us further suppose that 1+ M : Li, —k, for p=0 is smaller than all the other terms with p70. 


Furthermore, we suppose that the outer time scale of the neutron flux 7) = é/ (1 +M? ui - ky) 


is an order of magnitude greater than the time constant of evolution of variable w and the 
characteristic time of precursor decay. In this case, applying the second principle to link time 
scales, we can consider that both fields w of the feedback variable, and @ of the delayed 
neutrons. are in equilibrium with the neutron flux. 


With this assumption, the system of equations [D.1.4 a] reduces to the following: 

d 2,2 7 2y 

a (14+M?y2—k,)A, +k,K Sima, ~k,K Sn Ay +5S,,(t) 
mn m,n,q 


[D.1.5] 


Let us take as our time origin the instant in which all the external forcing terms Sp vanish. 
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Then, as the linear terms have negative coefficients, for a perturbation near enough to zero, the 


state of the system given by the mode amplitudes A, (t) will return to the origin of the space of 


mode amplitudes. 
Now, even when we assumed that 1+ a res = Ke is near zero, the coefficients of the linear 
parts of the equations for the other mode amplitudes are not near zero. This is always fulfilled in 


a bare core of standard dimensions. 
Let us suppose also that the external source excites mainly the zero-order mode amplitude, 


Ay (t), because S, (t) is much greater than S , (t) for p #0. Then the amplitude A, (t) will be 


dominant, and it is possible to uncouple it from the others mode amplitudes. 
We obtain the evolution equation for the uncoupled dominant mode (This applies for t<0O if the 
time origin is taken when the perturbation is already completed) 


d 
ae =-(1+ M743 ~ Ky JAy +k KT oq) Ap Hk KT ogo Ao +S,(t) Mons >0 5 Loo >0) 


[D.1.6] 
The others mode amplitudes evolve approximately this way, slaved by the dominant mode: 


<A, - —(14+ Mu -k)A, +k KT Ag —k KT o4p +S, (t) p #0 


[D.1.7] 
The others nonlinear terms are negligible. (The justification of this kind of assumptions can be 
found in Eckhaus, 1965) 
Now, the external sources are withdrawn in the instant t=0, and we have the initial conditions 


A,(0),A,(0),... with |A,(0)| much greater than the other|A, (0)| : 
Let us put S ,(t) =0. 


The resulting homogeneous equation: 


d 
apo —(14+M7u, —k )Ay+k, KT Ap — =k ,K* Toop 


k; Ay posi 


Has at least one real root A, = 0. It is the only real root if 1+ M7; —k, >. 
4k Tq00 


It is always stable. If it is the only real root, lim A, (t) =O and the higher mode amplitudes 
t-+400 


will be forced to approach to zero, as can be seen from equation [D.1.7]. 
272 
KT o00_ 000 
4k Ty 09 


Ao, 18 unstable and Ags is stable. The unstable mode amplitude is given by the formula: 


But if 1+ M75 —k, <1 there are two other real roots, Aoy andAgs. 


K Loo — ft L000 ~ 44 2! coo (1+ M7445 — ky ) 


Ag, = [D.1.8] 
2K k, Tay 


IfA (0) <A,4,,, then Ao(t) will return to zero dragging the other mode amplitudes to zero. 


But if Ay (0) >Ay,, A at) will grow tending to A, and dragging the other mode amplitudes 


to follow it. 
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We suppose that the time scales 7, = £ / (1 +M?* Lt, — k,) of the higher order mode amplitudes 
(p40) are an order of magnitude smaller than the time scale 7, of the dominant mode 


amplitude Ag(t). Then, after a short transient relative toz,, the slaved mode amplitudes will 


behave approximately as: 
A, (t)* 


So, the whole neutron flux will leave the set of states attracted by the zero-flux condition and 


KK Tj) Ag (t) KOT eg Ay lt) 


(ear =F [D.1.9] 


will tend to another stable steady-state flux. 
If the negative feedback is weak enough, the evolution equation for the dominant mode 


simplifies to: (A, = (1 + Mu, ~ky )Ag +k KT yoAg [D.1.10] 
t 


Besics 1+M?u 7 ky 
The threshold amplitude is given by A = k KI 
1 000 


The stable solution Agg now has moved away to infinity. Given an initial condition A (0), 


the solution of Equation [D.1.10] is given by: 


[D.1.11] 


Here oy = KT ony/ and @ = (14M ik Y/ +0 that = = Asa. 


Then, if Ag(0)< AouA o(t) approaches zero, and for ¢ large enough, it approaches zero as 


the damped exponential e 
But if A o(0) PA ig A o(t) will run away and for a finite time it will approach infinity. 
If we take ide as a parameter to vary, it is inversely proportional to the square of a characteristic 
linear dimension of the reactors core. 
k ce ae 
4k 1 00 


Then the zero flux will be the only steady state solution, and it will attract all the other states. 


If the core is small enough, su) will be big enough so that 1+M7*u5 -k, > 


So, no matter the amplitude of a perturbation, it will never cause instability. 
But if the core of the nuclear reactor is big enough to reverse the above inequality, a 
threshold of amplitude appears and with it, instability to finite perturbations is produced, 


even if the zero flux remains locally stable. 


The method allows us to calculate the global bifurcation of two branches, one corresponding to 


A 0, and the other to Ay, s» When ie reaches a critical value (Figure 2). 
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Fig. 2: Global bifurcation diagram for the equilibrium amplitudes of the dominant mode in terms 
of a bifurcation parameter inversely proportional to the square of a characteristic linear dimension 
of the reactors core. 


Now, as a dynamic example of possible oscillatory instabilities, let us consider nonlinear 
mechanical vibrations in a stretched shaft of uniform cross section. 


(D.2)-An example of dynamic instability 

Here we pose first the equations for transverse vibrations of a simply supported and non- 
rotating beam with a mid-plane stretching non-linearity. 

The vibrations occur always in one and the same plane, so the transverse deflection field is 
always parallel to the same axis orthogonal to the undisturbed beam axis. 

Then we consider the equations for transverse vibrations of a simply supported and rotating 
beam, also with a non-linearity due to mid-plane stretching. But now we need two transverse 
deflection fields, orthogonal to beam axis, that appear coupled in the equations of motion. 
Mid-plane stretching is in general a second order effect, so it is often neglected. In our case, it is 


a convenient nonlinearity to show the application of NMA methods. 


(D.2.1) Flexural vibrations without shaft rotation 
Let us consider first the transverse vibrations of a simply supported and non-rotating beam. Its 


length (when neither pre-compressed nor pre-tensioned) is/. It has uniform density p and 


uniform circular cross section of area A . Its Young modulus is F and its second moment of area 
is represented by/ . It is taken about an axis passing contained in the cross section and passing 
through its center, and sometimes called centroidal moment of inertia. There is internal 


dissipation of elastic energy: the normal stress o and the normal strain € are related by the 


equation: o=8[esr oe) [D.2.1] 
t 
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The beam movement is referred to an inertial system of Cartesian orthogonal 


coordinates x, y, z. The axis z coincides with the beam axis when the beam is at its equilibrium 
position: 0 < z</ 
The vibration occurs in the y, z plane, so the x direction is not involved, as it certainly is if the 


beam rotates about its axis z and also vibrates. 


The transverse deflection of the beam cross section center, measured from its equilibrium 
position, will be represented by w(t, Zi We will consider the possibility of having axial forces 
P with associated displacements v(t, z) transverse to the local beam cross section. These axial 
forces appear due to pre-compression, pre-tension or simply due to the lengthening of the beam 
when flexed. Besides, we assume a distributed external force (load) of linear density f (t, z) ; 


The following Raleigh’s equation (Meirovitch 1967; Meirovitch, 1997; Meirovitch, 2001) 
includes all the above-mentioned phenomena in a unified mathematical model: 


fake et 6° ot fake 
paZw(o)+E Hf ae w(t,z)+T, ams) Pl ae ap w(t,z) P ae w(t,z) =f (t,z) 


[D.2.2] 


Now, let us relate the axial force P with the displacements v ie z) and w(t, z) : 


An element of beam that when in equilibrium is in the interval [z, ztd z| and has lengthd z, 


after transverse displacement and elongation has length 
ds= (dz +dv) +(dw) 


Then, the longitudinal deformation A can be estimated by: 


2 2 2 
fA OE ray So af ol cape pe ee [D.2.3] 
dz Oz Oz Oz =2\ ez 
Assuming a linear relation between force and deformation P= EAA from [D.3] it follows: 
2 
P od 
Lye ae iy [D.2.4] 
Oz EA 2h.62 


Let us suppose now that the beam in equilibrium, after a longitudinal deformation related with 
the force P , has length] (1 + B) ‘ 


When there is pre-tension # > 0, when there is pre-compression # < (and if the shaft ends are 


fixed without neither pretension nor pre-compression, 2 = 0 


Then, if v(t,0) =0 and v(t,/) = 1, integrating [D.4] between z =O and z=/ we obtain: 
By Ie: \ 
pl=—i- } w| dz 
EA: : 25\ 0Z 


tire.) 
Sy P=EA d D.2.5 
) [a i(S / [ ] 
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From [D.2] and [D.5] we obtain the following nonlinear version of Raleigh’s partial differential 
equation: 


2 
2 


fa) o* 3° a 
aye Wz) + ee ( Smteden anes) 1G Bag V(b?) 


0z 
[D.2.6] 
ee ar aes en Oe 
—Cy pz iZ) dz age Witz) = P(t) 


I 
Here, by definition: =c,is the extensional wave velocity, ri =r,1is the radius of 


pA 


Ss Dm 


gyration and p(t,z)= is the external load per unit mass of shaft. The nonlinearity 


that appears in [D.6] results from mid-plane stretching and often it is a second order effect. 


Introducing the ansatz w(t, z) = >A, (t) sin( 5) in terms of the Eigen functions of the free 
n=1 
jRZ 


beam problem, in [D.6] and projecting onto the Eigen function sin 222) we obtain, after 
l 


several rearrangements: 


d dt Vek mj = 
| 2 
B ua +1 aj ey) ans 
- - P ale E’G| 74 
@,=C a [D.2.8] 26 0 =; ; [D.2.9] 
| e | Z 
1 a 
J p(tz)sin{ Ja 
P; (t) = v —<s [D.2.10] 
val ® J 
G P 
2 ee 
E'G [4 10; 
When B=0, w= oP and 2¢,@,=T,@; 80 ¢,= 7 
I+8[ 5 
When the first mode amplitude A io) is dominant, [D.7] reduces to a forced Duffing equation: 
a 
2 
a? d a 5 
2 3 
ap Milt) + 26101 Alt) + @7 A(t)}+ —~—4y Ail) = Pil?) [D.2.11] 


Considered as a one degree of freedom nonlinear oscillator, it is a hard oscillator because the 
slope of the elastic characteristic increases with the amplitude of oscillation. The properties of 
the solutions of this kind of damped and forced Duffing equation are detailed in Minorsky 
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(1983) from an engineer’s point of view and in Atlee Jackson (1989) from a mathematician’s 


point of view, amongst other places (additional references can be found in Zwillinger, 1997). 
4 


o{ a 
Cr - 
Defining q o) =——~_ ~~ and considering the harmonic external force of 


2 
l+r3 (=) 
frequency @ , p (t) = P19 COS (at) equation [D.2.11] can be written: 
d 2 
dt 


The external harmonic force has the following threshold amplitude: 


A(t) +260, A(t)+0} (A,(1)+aA3(1)) = p,9c0s(or) [D.2.12] 


(p jae a 
cae Ba 


When p,, is less than a threshold value ( Pio) , the amplitude A 


u 


[D.2.13] 


=A, of the steady 


1,max 
oscillation of the first mode as a function of the external frequency behaves as shown in the 
upper part of Figure 3. It differs from the purely linear resonance response only in the 
asymmetry of the curve. 


But if p,,> ( Bie) the resonance curve has three branches, two stable ones (branch ABC and 


branch DEF, represented by continuous lines) and one unstable (branch CD, represented by a 
dotted line). Only the stable branches correspond to actual oscillations of the first mode of shaft 


vibration. 


Figure 3: The ordinate coordinate axis represents the amplitude A, .,...= A, of steady oscillation. 
The abscissa coordinate axis represents the external frequency @. Above: p,,< ( Pip) Below: 


Pig> (Pr ), This qualitative figure corresponds to the resonance behavior of the nonlinear 


oscillator [D.12] for small damping. 
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If the external frequency increases slowly so that the amplitude A,increases following the 


branch ABC, at C the amplitude suddenly falls to the value corresponding to E (presents a 
decreasing jump). 

If the external frequency decreases slowly so that the response follows the curve FED, at D the 
amplitude suddenly increases to the value corresponding to B (presents an increasing jump). 


So, if p,.> ( Pi) a hysteresis phenomenon appears in the steady resonance response of the 


oscillator when the external frequency varies slowly enough: abrupt changes in the oscillation 


amplitude may occur. 


For small damping, besides the resonance when the frequency of the external force is near the 
natural frequency of the linear approximation, in principle resonance can appear also when the 
external frequency is near one third of the natural frequency (super-harmonic) and near three 
times the natural frequency (sub-harmonic). 

In principle also, if the external force oscillates with two or more frequencies, simultaneous 


resonances or a combination of resonance frequencies can occur. 


If the dominant amplitude A ,(r) is assumed to be known, the higher order mode amplitudes 
A, (t) (j =2,3,...) verify (approximately) a forced linear oscillator equation with an 


i mj? 
Cr [3 
A, 


undamped frequency Oo; + ve ("F] 1 (1 ) that varies with A(t) : 


l+r 7? 
rr +2 
a? d af jf 
2 2 
ae (1)+2¢, 0,7 A, (t)+ Os ay (1) A, (t)=p,(¢) [D.2.14] 
| ) 


When A,(t) is periodic, the solutions of equation [D.14] could show behaviors closely related 


with what is known as parametric resonance. 


When the dominant mode amplitude oscillates harmonically and the forcing term in the right- 
hand side of the equation vanishes, [D.2.14] is a Mathieu equation with linear damping. 
A detailed analysis of Mathieu equation can be found in Minorsky (1983) and Atlee Jackson 


(1989) amongst other places (additional references can be found in Zwillinger, 1997). 
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If the damping term in the Duffing equation [D.2.11] for the first mode amplitude can be 
neglected and the external forcing term vanishes, the angular frequency of the resultant 


nonlinear oscillation with initial conditions A,(0)=A, and <A,(0)=0 can be 
t 


approximately estimated by the formula: 


According to formula [D.1.15], the angular frequency increases with the amplitude of 


oscillation A, as always happens with hard oscillators. 


[D.2.16] 


This last formula allows a fast and easy analytical estimation of the period of the oscillations of 
amplitude A, corresponding to the solutions of the conservative nonlinear equation 


2 
“+03 f (x)=0 when f(x)=—f(-x) and f(x)=x+c,x° +c,x° +c,x’ +... (Detlaf and 
t 


Yavorsky, 1975) ° 


(D.2.2) Flexural vibrations with shaft rotation 
Let us consider an axially symmetric shaft whose lateral movement can be decoupled from its 


axial and torsional movements. 


> The well-known exact formula for the period is (A, pif ) ee (Jordan and Smith, 


} f(u)du 

1979) is in general cumbersome to be applied in analytical approaches, so during the past siécle a lot of 
effort was done to obtain more friendly approximate formulae (see for example: Mittal, 1978; Bravo 
Yuste and Martin Sanchez, 1989). 
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The spin velocity of this rotor will be represented byQ.. 
Now, unlike what happens when the axis vibrates but does not rotate, it is not possible to un- 


couple the bending motion in two independent movements in the principal planes y, z and x, z 


(Figure 4a and Figure 4b). 
\Z 


Figure 4 a (left, modified from Genta, 2005): In general, a rotating and deflected shaft is 
characterized by two angular velocities: its spinning velocity (2 around its curved axis and its 
whirl velocity around the Z axis. 

Figure 4 b (right, modified from Den Hartog, 1985): A cross section of the rotating and deflected 


shaft is shown. The intersection of the curved axis of the deflected shaft with one of its cross 
sections is S and the intersection of the original, non-deflected axis, with this cross section is B . In 
this case, the spin of the shaft around its center Sis suggested by a curved arrow in the upper left 
quarter of the figure. The whirling motion of S around B is shown as a broken circular curve 


with an arrow to suggest the sense of rotation, also in the considered case. 


Let us consider a rotating orthogonal coordinates system in the x, yplane of Figure 4 a, and 
with the same origin. The rotation velocity will be selected to be equal to the spin speed Q of 


the shaft. Now the lateral displacements u(t iz) and u,(t 4 of S will be referred to the 


rotating system. If w(t, z) and w v(t z) are its components relative to the fixed axes, we have: 


ry) fein satay eta 


w,(t,z)} | sin(Qr)  cos(Qr) || w,(t,z) 


The correction term that takes into account the effect of mid-plane stretching now depends of 


[D.2.17] 


the integral: 


f 2 2 1 2 2 
: | oy. | +) w, | Jaze i} Os | ee, de [D.2.18] 
21+, |\ Oz on 214 |\\ Oz oz 


(The equality follows from [D.2.17]) 
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From an application of the variation (energy) methods of continuous mechanics, it follows that 


the lateral displacements u At z) and u, (% z) verify the approximate coupled equations: 


‘tis: Xt ‘ : [D.2.19 al] 
—Cy —| [2m] +f «| dz sU,=P, 
21s OZ 0 ) 
or o* e° 
apt cif Sou sort | Ou 420 u, 
[D.2.19 b] 


Lue ee a? 
~r al [2] [Su] dz | ata Pa 


0 


Neither Raleigh’s rotation correction nor Timoshenko’s shear correction are considered here 


(Timoshenko, Young and Weaver 1974). The meanings of G ( extensional wave velocity) and 
ry (radius of gyration) are the same already explained in (D.2.1). 


When the spin Q2=0 the coupling between the lateral displacements fields u(t, z) and 


2 2 
1 | a é 

Uu, (t,z) disappear, in spite of the terms proportional to — | —u, | +| —u, | pdz, and 
21+ \\ Oz Oz 


the motion is produced in a plane that includes the z axis. Now, without rotation, the axes are 
X,Y, Z. 


If the axes x, y are rotated so that in the new position y, z coincides with the plane of motion 
1tl(a \ fa y 1¢(a ) 

—|5| —4, | +} —4, | pdz=—]| | —u, | dz and u(t, z) =0(0), so only one equation, 
ZN OZ OZ abe OZ 


that of w,(t,z)is left. 


Let us introduce the mode decomposition of the lateral displacement fields, corresponding to a 


simply supported beam: 


ty (t.2)= Ye Ara ()sin{ 2) [D.2.20 a] 

n=l 

u,(t.z)=>°A,, (1) sin( "5 ) [D.2.20 b] 
n=l 


From equations [D.2.19] and [D.2.20] we derive the following coupled equations for the mode 


4 +4 
' a 
amplitudes, with o; = e ie a : 
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2 


d d 
A, ; (t)+o; [a +T, £4, }ro%ay rr 


dt™ 
pants = [D.2.21 a] 
(= J {Sea (DA! fa, (1)= p,, (0) 
d* d d 
ae A,, (t) +0; [4,, + Ty ay A,; + Ay: +20 ae A,, 
[D.2.21 b] 


14 TE ][S mace 29}, W=0.00 


ofa fee ay' (uve ay S| 2 fess ay (Cea) 


ov 
a [D.2.22] 


dt*| A); dt| A,, 
2 +00 +00 A.. 
{SE} Sora. (eSemra, cco} 4 |=] Pe 
I n=l n=l 2j Pj 
Sa : ._|PA 9 
In [D.2.22] we have, with the diagonal inertia matrix M ( ji) = 0 A : 
fe 


My" (Duyeata))=| "2" sean al 


MYUNG Oo o]{p ol 


In this case the circulatory matrix N ( J 3Q) is identically zero. 


The linearized and homogeneous version of the coupled mode equations are the following: 


d* Bac aoe d 

Se Ai lt) +@jFa Ais +(@; -9 JA, ~225- Aa, =0 [D.2.23 a] 
2 

aA, j(t)+o77, <A, +(a@;-2°)A, ,+ 20 A,, =0 [D.2.23 b] 


Let us consider the asymptotic stability of the zero solution A, ; =0 , A, , =. This is the only 


rest point of the coupled system if @, # Q. 
This can be done (Jordan and Smith, 1979; Saaty and Bram, 1964) substituting in [D.2.23 a] and 
[D.2.23 b] the ansatz: 
A, ,(t)=A,,,(0) exp[ st] A, ,(t) =A, , (0) exp[s¢] 
These last expressions are solutions if and only if: 


P,(s;Q,7,)=s* +a,8° +a,8° +a,s+a, =0 [D.2.24] 
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2 
2 2 2 42 2 2 2 2 2 
a,=20;T, a,=2(a; a) )+ojr; a,=20;7,(0; =) ) ay = (0; a) ) 


According to the Routh-Hurwitz condition (Saaty and Bram, 1964; Attlee-Jackson, 1991) a 
necessary and sufficient conditions that assures that all the roots of [D.2.24] have negative real 
parts is the following: 


a,= 207, >0 [D.2.25 a] 
a, a, —a, =20; 7, [ +057; -? | >0 [D.2.25 b] 
a, , A, -(a; +p a,’) =40' | @; -? | >0 [D.2.25 c] 


The restriction [D.2.25 a] is always fulfilled when 7, > 0. The second restriction [D.2.25 b] is 


fulfilled if the third one [D.2.25 c] is verified. So, the necessary and sufficient condition for 


asymptotic stability of the rest point A, ; =0 ,A, , =0 of the linearized and unforced system is 
oO, > Q. In this case the coefficients of the polynomial that appears in [D.2.24] are all positive. 
But if @, < © there is no asymptotic stability of the rest point of the linearized system. 

When @, =Q, P,(s; Q= ;.T,) = (s° +20;7,5+0;T;)s" = (s+a?r,) s° 

So, equation P, (s; Q=0,, t,) =0 has two double roots: s=0 and s= -0; T, The stability 
of the linearized system when @, = Q is marginal. 

When 0@,<Q , P, (s; Q,7,) =s'+4,8°+a,s°+a,S+a has at least one negative 
coefficient, a,= 20; t; ( _ Q’) and if the difference between (2 and @ F is enough, 
a,= 2( a; —Q? ) + Q; T, will be negative also. 

Let us suppose that the spin speed Q is varied near @,. When @, = Q we know that s =Ois a 


double root. Now, we study the behavior of the roots that coalesce to this double root when 


Qtends to @ ), In this case, approximately: 


OQ? ~@ 2(Q?- 
pe : i) 1+ ws i [D.2.26 a] 
Oi Ty Oi Ty 
Y -a; 2(Q° —@ 
ak : uF iz i [D.2.26 b] 
Ot; OT, 
oa -—Q 2( @? —Q? 
If a, >, beingi = V-1 we have 5S, = ( — ) 1+i ( ; ) and S,is its 
Oi Ty O77 4 


complex conjugate. Both roots have a common negative real part. 


If a < (2 we see from equations [D.2.26] that both roots are real and (if © is near enough to 


@, ) both are positive. So, the linearized system is unstable and the mode amplitudes, once 
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crossed the frontier of marginal stability corresponding to the zero double-root, tend to drift 
exponentially without oscillating. 

When @, > Q. we know that the homogeneous nonlinear system is locally asymptotically stable, 
but the stability assessment when @, <Q. requires to study the full homogeneous nonlinear 
system: 


d 2 d 
— A, ,(t)+@; [a +t, £4,,}+0 A, ,—2Q— Ay; + 


[D.2.27 a] 


[D.2.27 b] 


The stability behavior of this homogeneous system and its forced version, with focus in the 
dominance of the lower order modes, will be studied elsewhere. 


(E)- Discussion and Conclusions 


In the example of hybrid modeling of power system given here, in section (C), both the electric 
generator and the turbine were considered from a lumped parameters point of view. Only the 
uniform shaft that connects the rotors of these two machines is considered from the point of 
view of distributed parameters. Bearings, fluid seals and their rotor-dynamics effects (including 
possible instabilities) must be integrated with the lumped mathematical models of the turbine 
and the generator. 

Besides, to simplify the mathematical model of the shaft, only torsional oscillations with linear 
dissipation were modeled. 

However, shafts are in general stiffer for torsion oscillations than for flexural oscillations. 
Lateral vibrations due to mass eccentricities, weak rubbing between rotating and static parts 
with hot spot production, friction induced lateral vibrations in synchronous generators, and 
many other causes often generate instabilities and unwanted vibrations in shafts assemblies of 
power plants (Stodola, 1924; Den Hartog, 1985; Lee, 1993; Genta, 1996; Machowski, Bialek 
and Bumby, 1998, Adams, 2000). 


The examples of lateral vibrations considered in sections (D.2.1) and (D.2.2) are highly 
simplified mathematical models of nonlinear mechanical vibrations, without and with rotation, 
respectively, in a mid-plane stretched shaft of uniform cross section with bearings without 
friction and with infinite stiffness. 

The representative beam can be considered as simply supported at both ends. 

The objective was to suggest how to apply NMA methods in suitable contexts, not to study a 
realistic model of a complex turbo-generator shaft. 


For certain purposes, a modal analysis based in a distributed parameter mathematical model of a 
shaft in a power system is not practicable at all. 
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Figure 5 shows a sketch of shaft which is an assembly of in series steam turbine rotors and an 
electric generator rotor, with the corresponding connecting rotors. 


Figure 5 (Modified from Adams, 2000): In the left, the high and intermediate pressure turbine 
rotors, in the middle the low-pressure turbine rotor and in the right the generator rotor. Turbine 
blades don’t appear in the draw. 


When flexible bearings and their lubrication must be considered in the model of a realistic 
mathematical model of a complex turbo-generator shaft, if the unbalance and the displacements 
are small enough, a linearization of the dynamic equations can be done, as it is done in the case 
of vibrating structures. 

In the linear approximation a system of differential equations with many degrees of freedom 
(perhaps with several hundreds of state variables) is usually posed: 


as d. . = 
M —74(t)+(D+G)— a(t) +(K+N)a(t)= P(t) 
The inertia M , damping D, gyroG , conservative K and non-conservative WN stiffness matrices 
are square matrices of order n (Adams, 2000; Genta, 2005; Lee, 1993; Meirovitch, 1967) 
In case there are no rotations in the system, the skew-symmetric gyro-matrix G vanishes, as 


well as the skew-symmetric circulatory matrix N and linear equations identical to those used in 
the description of structural vibrations are obtained. 


Now, let us return to NMA. When applied to the estimation of instability thresholds (the sizes of 
a perturbation from a steady state that causes instability) the NMA methods can be considered a 
nonlinear extension of Liapunov’s first method. 


In the case of restricted NMA the equations for mode amplitudes A ; appear uncoupled in the 


linear approximation. For example: 


SA, = AAD TA At lini Ae A Ay to (jG 108s 
k=1 1=1 k=1 1l=1 m=1 


So, in the linear approximation A j (t) behaves as exp [4 j r| : 


When mechanical vibrations in continuous systems are involved is also possible for the 
uncoupled modes in the linear approximation to behave as linear oscillator, sometimes with 


damping and forced by the projection p, (1) of distributed external loads onto mode j 


2 
Fade J A) +26,0, 54, +03 4,=p,(0) 

AS was mentioned in section (D.2.1), this linear approximation appears when modal analysis is 
applied to a nonlinear Raleigh model of a transversely vibrating uniform beam with internal 
dissipation and nonlinearity resulting from mid-plane stretching. 

In (D.2.1) we saw the possibility of a hysteresis phenomenon related with adiabatic changes in 
the frequency of a harmonic external force driving the first mode uncoupled (Duffing forced 


oscillator) when the force amplitude is greater than a certain threshold. 
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However, formula [D.1.13] shows that the threshold amplitude increases when the strength of 
the nonlinearity @ decreases. So, when the nonlinearity is a mild one, like that due to mid- 
plane stretching, we can expect that the amplitude of the external force does not exceed the 
threshold amplitude. 


In general, when modal analysis is applied to continuous mathematical models of vibration in 
rotating systems, mode amplitudes appear coupled in the linear approximation. 
For example, in a continuous mathematical model of a uniform rotating and vibrating shaft, for 


each mode j ( j =1,2,...) two mode amplitudes Ai, and A); appear linearly coupled, as was 


seen in (D.2.2). The corresponding differential equations are of the following general form: 


me ba +(D(J)+ aaa) Slt H(K (i) Tee) WAR 


Here the matrices are in general dependent of the considered mode. Besides, the elements of at 


least two these matrices (G( jQ) and N ( 7Q)) are functions of the spin velocity Q of the 
shaft. 
my, (i) i 


ma, (j) My'( 7) 


bas 82 (42) 


. ; &y = d,,( ') dy ( ) 
mani D(J}+0(20)=| 0) 42())"sn:2) 8cC2) 


symmetric positive definite damping matrix D( ji) with a skew-symmetrical gyro-matrix 


The inertia matrix is symmetric and positive definite: M ( 7, ) = 


is the sum of a 


G( J 3Q) (mainly due to gyroscopic effects in rotating systems). 


= ee oy | ku) Ra) | fn (32) m2 (432) | 
The matrix (ema) =| ir hae an He is the sum of a 


symmetric stiffness matrix (due to conservative forces in the linear approximation) with a skew- 


symmetric circulatory matrix (mainly due to non-conservative forces in the linear 
approximation). 


Pij\2) |. bat’ Fo 
The vector q Js a forcing term due to the projections of distributed external loads onto 


P2,(t) 


mode /. 


The analytical approach of NMA works when there are a few low order dominant modes that 
enslave the higher order ones as assumed in case of the order parameters defined in Synergetics. 
Besides the kind of analytical justification of dominance given by Eckhaus, this dominance can 
be checked using digital simulations done from the properly truncated equations of the coupled 
mode amplitudes (Sudrez-Antola, 1997). 

An important additional assumption, underlying this dominance of the first modes, is the 
possibility of ordering the real parts of the eigenvalues and that they result to be widely 
separated. 
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To simplify the following discussion, let us suppose that for certain values of system’s 


parameters the eigenvalues are all real, negative and can be ordered: O>24,>/1,>A, >... and 


that the evolution equation are the following: 


d +00 +00 +00 +00 +00 


rr =A, A, +>) > TA Art > jim de A Ay te (7 =1,2,3,...) 


k=1 1=1 k=1 1[=1 m=1 


If € gives the order of a perturbation from the rest state and 6; is a multiplier, we introduce 
new mode amplitudes a, (t) through the relation: A, (t) = E0, a; (t) 

Now, some parameter is modified and approximates to a certain critical value whered ,=0 . 
Then we assume that while / , tends to zero the relative ordering of the eigenvalues remains the 


same and they don’t tend to zero. 


A ; ; ; 
So, let us take ¢=—A, and 6 — ri so the evolution equations for the new mode amplitudes 
J 


y ie 

d +00 +00 

are: a =A,a,;-A, EEA Joa. # C7H123345) 

k=1 [=1 A, 

Near the point of marginal stability, where ¢ =—/ , tends to zero, the only non-vanishing term 
+00 +00 As A+ 

in DT age aa, is Ti, TA a,a,=I1,,,4,qa, and an analogous results with 
k=1 [=1 1?" 


the terms of ni order in the products of the new mode amplitudes. 


AS consequence, near marginal stability: “a, = A, a; a, —A, I @ +... (7 =1,2,3,...) 
; 4 . d 2 . 
Returning to the original mode amplitudes: aoe =A, A, +1, A, +... (J =1,2,3,...) 


Now the equations for the higher mode amplitudes j = 2,3,.... are coupled only to the first one, 


so their dynamics are forced by the first mode amplitude, which is then the dominant one and 


d 


varies according to: ai —A, =4,A 4h, A +4, 


A; + 


If they exist, the threshold amplitudes appear in this last equation. 


For example, being 2, negative, if J,,, is positive and higher order terms, beginning by the 


4, 


third order, can be neglected, there is a threshold amplitude A 


lu 
111 


If the initial value of the dominant mode, A, (0) <A,,, thenA ale ) tends to zero and drags to 


lu’ 
zero the other mode amplitudes. If A, (0) >A,, then A, (t ) tends to grow and drags enslaved 


with it the other mode amplitudes. 


Although Eckhaus original approach to NMA, like the approach to order parameters in 
Synergetics, seek a solution valid near marginal stability (when the real part of dominant 
eigenvalues is near enough to zero) the example of the nuclear reactor static instability and the 
results of restricted NMA applied to the study of the threshold of electric excitation of 
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biological tissues (Sudrez-Antola, 1997) suggests that marginal stability as such is not a 
necessary precondition. 

As marginal stability is not needed, it was possible to find a static global bifurcation that 
predicts a possible power runaway in a nuclear reactor, given suitable circumstances. 

This result could be of some interest in relation with the physics of subcritical multiplying 
systems driven by particle accelerators. 

Another interesting case, always in the field of nuclear reactor dynamics is the hard excitation of 
the so-called xenon oscillations in a PWR (Sudrez-Antola, 2005 b). Here the threshold 
amplitude of the above-mentioned hard oscillations is estimated from the coupled nonlinear 
evolution equations of mode amplitudes and corresponds to an unstable limit cycle that 
surrounds a stable steady state. 


However, threshold estimation by NMA presents certain important limitations. 

The form of the initial perturbation should be dominated by the eigenvectors of the dominant 
modes. 

It could be possible that another kind of perturbation could cause instability, perhaps even at a 
smaller threshold amplitude. 

The success of NMA methods depends on what could be considered as a kind of continuity 
between the linear and the nonlinear domains of systems dynamics. But not all interesting 
behavior is under the scope of this kind of continuity. 


The other disadvantage is that if there are not a few dominant modes the analytical approach 
probably will not succeed, and a numerical calculation should be undertaken. 

In this point it is necessary to assess if it is better to leave NMA equations and use a computer 
code for an ab-initio discretization of field equations. 


In any case, restricted NMA methods that allow an estimation of the response to finite- 
amplitude disturbances using an extension of Liapunov’s First Method are straightforward in 
principle. 

They can provide a considerable amount of information, of value by itself and for the design of 
digital simulations applying computer codes corresponding to more complex and realistic 
mathematical models of the studied system. 


The calculations needed to apply the asymptotic methods of NMA may be long and tedious. 
However, powerful symbol manipulation packages are now available to develop complex 
symbolic calculations in the computer. 

Of course, these tools were unavailable when the original approach to restricted NMA was 
constructed. As consequence, daunting pencil and paper calculations had to be undertaken 
during the sixties and the seventies to apply these methods. Then, it was natural that people 
shifted to other analytical tools easier to use, and to numerical methods, applied from the very 
beginning. ° 


© Nevertheless, both new and useful results can be obtained making digital simulations with the ordinary 
differential equations obtained from partial differential equations applying NMA. A famous example is 
the numerical discovery of chaos, at the beginning of the sixties, during numerical simulations of the 
dynamics corresponding to Lorentz equations. These nonlinear ordinary differential equations are the 
equations of evolution of the coupled mode amplitudes obtained from a very simplified model of an 
atmosphere (Attlee-Jackson, 1991; Nicolis, 1995). 
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The situation now is different. Even if present day symbolic packages may be still crude for 
certain specific manipulations, the processes of improving them are going on very fast. 
Furthermore, the contemporary tools of nonlinear science and a very significant amount of 
experience in applying them is already available. 


(F) Appendices 
(F.A)- Weighted Residuals and Galerkin Method 


The methods of weighted residuals are used to convert a partial differential equation with time t 
and space 7 as independent variables into a system of ordinary in the time variable only. 


Let u (t, ia ) be a field that is a solution of the partial differential equation: 
O| u(t.) ]=f (t,7) [F.A.1] 
Here O is a given operator, dependent of time and space variables, and f (t,7) is a given 


source term. Both the operator and the source term are defined in a spatial domainD 
for 0 < t < +00. Suitable initial and boundary conditions are established to have one and only 
one solution for the partial differential equation. 


Let u (t, r ) be an approximate solution. Then, its residual is defined to be: 


R(t,7) =O[a(t,7) |-f (47) [F.A.2] 
We could think that a good approximate solution is one that makes the residual small. If the 


residual is zero in every place in the domain D and for every time instant, then w (t, r ) would 


be the solution of the problem. But, in general it would not be the exact solution we are 
seeking for, so we must define what we mean with “the residual is small”. 


Now, let us consider a denumerable infinite and complete set of space functions \d, (F )} in the 
space of functions to which the solution u(t ,x) belongs for every fixed time instant. These 
?, (F ) are known as base functions and are chosen so that they satisfy the boundary 


conditions of the partial differential equation problem. However, they need not be orthogonal. 
The set of base functions being complete, we know the solution can be expressed as a series: 


u(nF)=DA (t)-¢, (7) [F.A.3] 


Usually we can’t find the whole infinite series, only a finite approximation, and so only an 
approximation to the full solution: 


i(t,7) =D) A, (t)-¢, (7) [F.A.4] 


N 
k=l 


The Bubnov-Galerkin method (or simply the Galerkin method, as is often referred to) is a 
specific procedure to minimize the residual obtained by substitution of [A4] in [A2]. 


The projections of the residual onto each of the functions ¢, (F) are restrained to be zero: 


[RO7)-¢(7)-dF¥=0 kK =12,...N [F.A.5] 
D 


We obtain thus N equations in the N unknown functions A, (t), A, eee Ay (t). 
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It the operator O implies taking time derivatives of the field (t, r ) , then due to the assumption 


[A4], this means that these time derivatives are taken on the’ time 


functions A, (t), A, (egAy (t). As consequence of the integration over the space variables 


done in the restriction [4], these space variables disappear and only the functions 


A, (), A, (t),..-, Ay (t) and their time derivatives remain: a system of ordinary differential 


equations is thus obtained. 
The initial values of A, (t), AS 3 Remy (t) and, if necessary, of their low order derivatives, 


stems from the initial conditions that were established for the original field u(t, X) whose 


approximations are being constructed. 
Example 
Let us consider the following nonlinear wave equation in one spatial dimension: 


oe o 0e : 
F aless)=e Zu(eye[ SZ ule) [Ex 1] 


N 
We take D = (o,1] and substitute #(t,x) = » Ai (t)@, (x) in the residual equations [F.A.4]: 


Wf a? 18. ale. ; 
} azilx)-¢ arneors[aleeet) d,(x)dx =0 [Ex2] 


cy (Sa, (9) (44 oF Se (x) Sa, (x)¢, (| =0 


(1 15 wt [Ex3] 
l d’ d’ d’ 
Now, | # (x), (x) dx, joo (x )d, (x )dx and iS b,(x) sa 


for k,J =1,2,...,N , as well as the wave velocity c and the damping parameter ¢ . 


o, (x )d, (x ) dx are known 


[ Ai) 
A,(t) a 
So, introducing: the column vector a(t) =| , the square matrices of order N, A with 
LAy (+) 
Ly 
sieeicatel d, () (x)@(x)dx, and B of elements [5 (x)dx k,1=1,2,....N, and the 
0 » 


= 2 2 
third order tensor C with components rie 6) (Sa (x) (sa 


2 
NaS 


(j,k,/ =1,2,...,N ), we obtain this vector nonlinear ordinary differential equation: 
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2; 
AG a(t) 5 LalteCeSa—c! B-a()=6 [Ex4] 


In [E4] the symbol ¢ represents the scalar product (Brand, 1947). 
To go further we need to specify the boundary conditions in x =Qandx=/, and the initial 


conditions u (0, x) and Ky (0, x) : 
t 


(F.B)-A lumped parameter model of torsional oscillations 


Let us return to the simplified hybrid mathematical model of the exchange of mechanical and 
electrical energies built as a generalization of the classic Stodola (1917) vibration model in 
section (C) “Hybrid distributed parameters-lumped parameters mathematical description of 
power systems” of this paper. 


Our purpose now is to derive, from the hybrid model, a lumped parameters model able to 
include some details related with the interaction of the synchronous generator with power lines 
and loads, so that some aspects of sub-synchronous resonance phenomena can be studied. 

The following lumped parameters model also generalizes the above mentioned Stodola model 
but is easier to study than the hybrid model. 


; 20 
Let the maximum angular frequency involved in the dynamics be@,,,, =—— . 


min 


Then 


is the quotient between the shaft length and the length traveled by a perturbation 
Cr 


min 


that moves with speed c, during an interval of time of duration T,,,, . 


Now we define: é=27n [F.B.1] 


min Cr 
Introducing the dimensionless time variable 7=@,,,.¢ and space variable ¢ =F? equation 


[C.3] can be recast as follows: 


3 


Co 
2 
é 

Or 


a” O(t,g)+o@ 


(7.5)= aa me TBO 


(z,é) [F.B.2] 


If the dimensionless number277 = &€ <1, wave propagation can be neglected, then: 


Cr 


min 


2 3 


P 
ae) tas 


O(t,x)=0 [F.B.3] 


As consequence, being A (t) an (for the present) arbitrary time function: 


0 ole 
By Ot) + Fa Ban =A) [F.B.4] 
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Integrating equation [4] relative to the position along the shaft axis betweenx =O andx=L: 


{2 os+r, <o.sfar-(2 Orb), d é 0-80) 400 


Thus: M,(1.x)=G4, (2080). a [t0-4 oh 


dt 


[F.B.6] 


So, if both M,,, (t) and M,., (t) are assumed to be known, the following lumped parameters 


mathematical model with only two degrees of freedom is obtained: 


S) a5 (t)=-M,, (t)+GJ, (2240).. : é aaa) [F.B.7 a] 


’ dt dt 


a @, (t)=-GJ, (2240).. eral (t)  [F.B.7b] 


‘dP “dt L 
Let us define the twist angle: Ag=¢, (t) —@ (t) [F.B.8] 
, GJ,{ 1 1 
And the twist frequency: oO, = tS [F.B.9] 
L\oO, 9, 


From equations [7 a] [7 b], considering the definitions [8] and [9], we derive a twist equation: 


d* d 1 1 
a bolt) +T,O, oot es Ag= @ Men (1)+ 5 Mm (t) [F.B.10] 


g t 


During steady operation the prime-mover (turbine) torque M pm Temains constant: 


m 


M j,=™M ym o In turn, the electromagnetic torque on the generator rotor 7 M.,,, (t) has a 
constant component equal to the prime mover torque M,,,=M pie =M, and a small 


oscillating component whose frequency is twice the synchronic one @, : 
6M, (t)=6M,,,,cos(2@,t+¢) 


These two terms of the torque, M, and 6M, (t) appear due to the following reasons. 


7 Both the rotor and the stator of the electric generator are assumed to be smooth (Woodson and Melcher, 
1968). 
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The whole electromagnetic torque is proportional to the DC rotor current, to the AC synchronic 
harmonic oscillating current in each stator winding, and to the derivative of the mutual 
inductance between the rotor winding and each stator winding, relative to the rotation angle. 
This last derivative oscillates harmonically with the synchronic frequency. 

The product of the AC stator currents with the mentioned derivatives of the mutual inductance 
gives a constant (DC, non-oscillating) torque added to oscillating (AC, harmonic) torques of 
twice the synchronous frequency (Woodson and Melcher, 1968). 

During steady operation, the average value of this torque of electromagnetic origin, over a 
period of alternating current is zero. 

But its instantaneous value oscillates changing sign. 

Taking into consideration these remarks, we see that during steady operation equation [F.B.10] 
reduces to the following: 


d 1 1 1 
—_Ad(t)+7, 0 —Ad+a? Ad= + M,+—oM, ,cos(2@.t+ [F.B.11] 
o(t) d 4 dt d A d © 2) 0 2) em,0 ( Ss 2) 


g t t 
After a possible transient decay, the solution of this equation is: 

Ag=A¢ ot Og(t) 
Here Ag ,=——M, is constant twist and 6 g(t) oscillates harmonically with angular 


P 


frequency 2.@, and zero average value: 


5 9(t) -( se) : cos(2@,t+9-a(a,)) 


2, 
oO, —4@,) +7, 0, 


2 

T,0 
Here: tga(@,)=—-+~ 
QO, — 40; 


In general, mainly due to the rotational inertia of the common axis with its attached masses, the 


characteristic angular frequency @, is at least an order of magnitude smaller thang, . If the 
natural frequency of the damped system is near@, , the fluctuation of the electromagnetic 


moment due to the purely reactive components of the load is filtered by the mechanical system 
and almost does not influence its steady (regime) angular velocity. 

It causes small alternating mechanical stresses on the axis that are not negligible from the point 
of view of the strength of the materials and the mechanics of the fracture, because they can 
cause fatigue growth of defects. 

All this must be considered, first during mechanical design, and then during preventive 
maintenance of the machine, but this small fraction of mechanical energy dissipated during shaft 
deformations will be neglected here. 


However, as was said above in section (C), a short-circuit due to a system fault, or a sudden and 
significant load change in a compensated power line could produce an exchange of energy 
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between the shaft and the inductive and capacitive elements of the electric circuit, generating 
very low frequency mechanical torques on the generator rotor. 


Let us suppose (Machowski, Bialek and Bumby, 1998) that the compensated power line feeded 
by one phase of the synchronous generator, after the fault or load change is represented by an 
in-series combination of a linear inductance L, a linear capacitance C,, a linear resistance R 


and a source of voltage of amplitude U, (the induced voltage in the stator windings 


corresponding to the synchronous angular frequency q@, before the perturbation). The 


Ca 


AY 


impedance after the fault is thenZ = R+ i Lo,- | 


peneg = JE (assumed to be less than one),Q, = andQ= 1-67 Q.,, 


1 
Tie 
being Q<a@,. 

The phase current i(t) after a short circuit can be expressed as the sum of a steady synchronous 


oscillating term i, (t) proportional to the amplitude of the phase voltage source and an 


exponentially damped transient term i,, (t) oscillating with the angular frequencyQ : 


i(t) =i, (t)+i, (t) = gin t+@, ) +45 exp[-¢ Q., t]sin(Qr + 2,) [F.B.12] 

As said before, the whole electromagnetic torque is proportional to the DC rotor current, to the 
oscillating current in each stator winding (now the sum of the synchronous AC steady term and 
the transient term that appear in the formula [12] above), and to the derivative of the mutual 
inductance between the rotor winding and each stator winding, relative to the rotation angle. In 
a first approximation we can assume that this last derivative oscillates harmonically with the 
synchronic frequency. 

The product of the stator currents with the mentioned derivatives of the mutual inductance gives 
four oscillating terms added together to give the whole electromagnetic torque. 


-A constant (DC, non-oscillating) torque. 


-A torque 0M cos(2@, t+, ) oscillating harmonically with twice the synchronous 


em,0 
frequency. 


-Another torque 6M,,, «1.5 exp[-¢ Q, t]| sin [(2, 20 )t + 2, | with the sub-synchronous 
frequency @, -Q. 
-A last one OM,,, cys) €XP [-¢ O. t] sin (2, +Q )t+ Py | with the super-synchronous 


frequency @,+ Q. 


If these low frequencies (2 are such that @, —Q is near one of the (lowest) natural torsional 


frequencies of the shaft, strong resonant vibrations may appear in the shaft assembly (Figure 
F.B.1). 
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electromagnetic torque 


oscillatory response 
of the shaft assemblv 


“0 02 04 06 00 1 12 i4 16 Tr: 2 
time [sec] 


Fig.F.B.1: (Above) Transient sub-synchronous torque 6M .)., ¢, ., exp[-¢ oO t]sin [( a, -Q e] on 


the rotor generator, of electromagnetic origin, due to a system fault or a sudden load change. 
(Below) Resonant growing torsional oscillations induced by the transient torque on the rotor 
generator and energized from the prime mover (example modified from a Siemens report). 
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